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ABSTRACT 

Context. Colliding flows are a commonly used scenario for the formation of molecular clouds in numerical simulations. Due to the 
thermal instability of the warm neutral medium, turbulence is produced by cooling. 

Aims. We carry out a two-dimensional numerical study of colliding flows in order to test whether statistical properties inferred from 
adaptive mesh refinement (AMR) simulations are robust with respect to the applied refinement criteria. 

Methods. We compare probability density functions of various quantities as well as the clump statistics and fractal dimension of the 
density fields in AMR simulations to a static-grid simulation. The static grid with 2048 2 cells matches the resolution of the most 
refined subgrids in the AMR simulations. 

Results. The density statistics is reproduced fairly well by AMR. Refinement criteria based on the cooling time or the turbulence 
intensity appear to be superior to the standard technique of refinement by overdensity. Nevertheless, substantial differences in the flow 
structure become apparent. 

Conclusions. In general, it is difficult to separate numerical effects from genuine physical processes in AMR simulations. 
Key words, adaptive mesh refinement - hydrodynamics - turbulence - thermal instabilities 



1. Introduction 

Computational fluid dynamics in astrophysics relies on numeri- 
cal methods that are capable of covering a huge range of scales. 
Apart from smoothed particle hydrodynamics (Monaghan 
|1992| ), adaptive mesh refinement (AMR) has been applied to 
variety of problems. This method was developed by Berger & 
|Qliger| ( ^9g4] ) and |Berger & Colella| ( [l989| ). Among the widely 
used, publicly availa ble AMR codes for astrophy sical fluid dy- 
namic s are FLAS H ([Fryxell et al ||200Q| ), Enzo ( [O'Shea et aL 
2004 ) and Ramses (Teyssier 2002JT Althoug h there are compar- 
* O'Shea etal.2005; 



ative studies of AMR vs. SPH (for exam ple, 
|Agertz et"aLl|2007[ |Commercon et a l. 2008), the degree of re- 
liance of AMR in comparison to non-adaptive methods has met 
only little attention so far. 

Especially for turbulent flows, it is a non-trivial question 
whether the solutions obtained from AMR simulations agree 
with the correct solutions of the fluid dynamical equations at a 
given resolution level. For this reason, we systematically com- 
pare AMR and static-grid simulations for a particular test prob- 
lem in this article. We chose a scenario that has been inves- 
tigated in the context of molecular cloud formation, namely, 
the frontal collision of opposing flows of warm atomic hydro- 
gen at supe rsonic speed ([Heitsch et al.|2006[|Vazquez-Semadeni| 

et al.|2007l|Hennebelle & Audit|2007al IHennebelle et al.|2008l 2. Numerical methods and simulation setup 



multi-phase medium evolving in these simulations is highly 
resolution-dependent, and numerical convergence is seen only 
at resolutions well above 1000 2 . In three-dimensional simula- 
tions, such high resolutions are infeasible if static grids are 
used. C onsequently, Hennebelle et al.|p008| ) and |Banerjee et al.| 
(2008) applied refinement by fixed density thresholds and re- 
finement by Jeans mass, respectively, in their three-dimensional 
high-resolution AMR simulations. 

In this article, we consider two-dimensional colliding flows 
without self-gravity and magnetic fields for a systematic com- 
parison of AMR simulations to a reference simulation on a static 
grid. We analyze both statistical properties and the morphology 
of the gas fragmentation due to the cooling instability. This work 
is organized as follows: In Section 2 the numerical methods are 
described and the setup of the simulations will be presented in 
detail. In Section 3, we compare the results from the different 
simulations. Section 4 concludes this paper with a summary of 
the main results and general remarks on AMR. 



Walder & Folini 2000). Because of the cooling instability at den- 
sities ~ 1 cm" J and temperatures of a few thousand Kelvin, the 
gas becomes highly turbulent at the collision interface. Since 
the instabilities develop on length scales much smaller than the 
integral scale, this problem is computationally extremely de- 
manding. The two-dimensional resolution study of Hennebelle 
& Audit (2007a) showed that the properties of the turbulent 



The simulations presented in this article are accomplished us- 
ing the open source code Enzo ( [Bryan & Norman|1997[ [(TShea 
et al. 2004| ). The compressible Euler equations are solved by 
means of the staggered grid, finite difference method Zeus 
dStone & Norman|1992a|bl |Stone et al.|T992]). We included the 
cooling function L defined by Audit & H ennebelle (2005]) in 
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Fig. 1: Phase diagramm of log(P) vs. log(ft). The solid black 
curve indicates the equilibrium curve defined by L - 0, and the 
straight lines (orange in online version) are the isothermals for 
7000 K and 50 K . The contour shading indicates the probability 
density. 



these equations: 
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Fig. 2: Contour plot of the mass density after 5 Myrs of evolution 
in the static-grid simulation. The density scale is logarithmic. 



d 

dt 



-p + u • Vp = 
p + V(pu ®u + P) = 



dt 



e + V[u(pe + P)] = -£(p,T). 



(1) 

(2) 
(3) 



The primitive variables are the mass density p, the velocity u and 
the specific total energy e of the fluid. The total energy per unit 
mass is given by 



1 



u 2 + 



(4) 



(7 "DP 

where y is the adiabatic exponent and the pressure P is related to 
the mass density p and the temperature T via the ideal gas law: 



pk B T 
fim H ' 



(5) 



The constants ks, V and m# denote the Boltzmann constant, the 
mean molecular weight and the mass of the hydrogen atom, re- 
spectively. The gas is assumed to be a perfect gas with y = 5/3 
and fi = IAjtih- 



The cooling function of Audit & Hennebelle (2005) in- 
cludes the cooling by fine-structure lines of CII and OI, fur- 
ther the cooling by H (Lya-line) and the electron recombination 
onto positively charged grains. The heating is due to the photo- 
electric effect on small grains and polycyclic aromatic hydro- 
carbons (PAH) caused by the far-ultraviolet galactic background 
radiation. For more information about this cooling function see 
|Wolfire et aL] (^995) [2003] ); |Spitzer| ft978) ; |Bakes & Tielensl 
" 1994j > and |Habing| (l968| ). The pressure-equillibrium curve re- 



sulting from the cooling function is plotted as black curve in 
Figure [T] For the numerical solution of the fluid dynamical equa- 
tions, we used the radiative cooling routine implemented into 
Enzo. For each hydrodynamical time step, the state variables are 
iterated over several subcycles, and the resulting total energy in- 
crement for the whole time step is added. 



For our numerical study, the tw o-dimensional setup of |Audit 
|& Hennebeflel ( [2005] ) and |Hennebelle & Audit| P007b| ) was 
adopted with small modifications. The initial gas content cor- 
responds to the warm neutral material (WNM) in the ISM, i. e., 



the temperature is T = 7100 K, the pressure is P = 1 x 10" 13 
erg cm -3 and the number density of neutral hydrogen is n - 
0.71 cm -3 . From the left and the right boundaries, warm gas 
with identical thermodynamic properties flows into the compu- 
tational domain, where the cosine-shaped inflow velocity pro- 
file is modulated with small perturbations, realised by randomly 
shifted phases. These phase shifts are kept constant for the 
different simulations, so that initial conditions are exactly the 
same for all runs to ensure comparability. Following Hennebelle 
et al. ( |2008| ), the top and bottom boundary conditions are peri- 
odic. The physical dimensions of the computational domain are 
20 x 20 pc. The two inflows of gas collide in the middle of the 
domain. The supersonic collision causes a steep rise in the gas 
density that triggers the thermal instability, and gas undergoes a 
transition into the phase of the cold neutral material (CNM) in 
the ISM. In this phase, the gas has temperatures i n the range 30- 
100 K and number densities within 20-50 cm" 3 dFerriere|200i] . 
The thermal instability produces highly turbulent structures (see 
Figure [2]) with Mach numbers up to 20. The challenge for AMR 
is to track these turbulent structures as accurately as possible. 

A reference simulation was run with a static grid of 204 8 2 
cells. Then the same setup was evolved in AMR simulations with 
a root- grid resolution of 128 2 cells and 4 levels of refinement. 
The resolution between adjacent refinement levels increases by 
a factor of 2. Hence, the effective resolution at the highest level 
of refinement is 2048 2 . In these simulations, we employed three 
different types of refinement criteria: 

1. refinement by overdensity (OD), 

2. refinement by cooling time (CT), 

3. refinement by rate of compression and enstrophy (RCEN). 

The first two criteria are widely used in astrophysical AMR sim- 
ulations. For refinement by overdensity, the mass density must 
exceed the initial density on the root grid by a certain factor. This 
overdensity, in turn, defines the initial density for refinement at 
the first level of refinement and so on. We chose three differ- 
ent values for the overdensity factor, namley, twice the initial 
density (default OD), as well as three times (OD-3) and four- 
times (OD-4) the initial density. For criterion CT, on the other 
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Fig. 3: Pdfs of the number density n (upper panel) and the temperature T (lower panel) in log-log-scale for the different AMR-criteria 
(black curves) compared to the static grid simulation (red curves in the online version). 



hand, refinement is triggered for a grid cell if the cooling time 
f"cooi := P/(y- l)p|X| becomes less than the sound crossing time 
over the cell width. Refinement by the rate of compression and 
the ens trophy uses yet another technique. It was introduced by 
Sch midt et al.| ( |2009 ) for the simulation of supersonic isothermal 
turbulence. The control variables for refinement are the enstro- 
phy and the rate of compression. The enstrophy is given by one- 
half of the square of the vorticity, while the rate of compression 
is defined by the substantial time-derivative of the negative di- 
vergence of the velocity. The expression used by Sch midt et al. 
( 2009) to evaluate the rate of compression (see equation (12) 
in this paper) is easily generalized to the non-isothermal case, 
where the speed of sound is not a constant. To trigger refine- 
ment by RCEN, dynamic thresholds are calculated from statisti- 
cal moments of the control variables: A grid cell is flagged for 
refinement if the local fluctuation of a control variable becomes 
greater than the maximum of the average and the standard de- 
viation of the variable. On the root grid, averages and variances 
are computed globally, whereas averaging is constrained to indi- 
vidual grid patches at higher levels of refinement. 

For comparison of the simulation results, we calculated 
probability density functions (pdf) of several quantities. To ana- 
lyze the gas fragmentation in each simulation, we adapted the 
clumpfind algorithm implemented by Pa doan et al.| ( |2007| ) to 
non-isothermal problems. The algorithm identifies the small- 
est dense regions that fulfill the Jeans criterion for gravitation- 
ally unstable gas. Since the clump samples found on the two- 
dimensional grids used in our simulations are insufficient for the 
calculation of clump mass spectra, only the total number and the 
mean size of the clumps are used for quantitative comparisons. 
In addition, we computed the fractal dimension of gas at den- 
sities higher than n = 20 cm -3 (corresponding to the minimum 
density of gas in the cold pha se) by means of the box-counting 
method ( |Federrath et al.|2009| ). 



3. Results 

Due to the gradual accumulation of gas in the simulation do- 
main, no strict statistical equilibrium is approached. For this rea- 
son, we evolved the flow until noticeable small-scale structure 
has developed and the separation of the gas into two phases has 
emerged. As shown in Figure [T] two distinct phases are found 
at time t = 5 Myrs. At this time, the central flow region is 
in a turbulent state (a contour plot of the mass density of the 
gas is shown in Figure [2]). Thus, we carry out our analysis for 
t = 5 Myrs. While the main fraction of the gas is situated in the 
warm phase with temperatures between 5000 and 10000 K and 
low densities ~ 1 cm -3 , the cold gas with temperatures between 
30 K and 100 K and densities in the range 30 - 350 cm -3 can be 
found close to the equilibrium curve. 

The pdfs of the mass density and the temperature obtained 
from different AMR simulations are plotted in Figure [3] In prin- 
ciple, all refinement criteria reproduce the distributions found in 
the static-grid simulation quite well, although there is a trend of 
slightly more cold gas at the cost of warm gas. The discrepancy 
is more pronounced for refinement by over-density (OD) com- 
pared to the other criteria, and it becomes worse for the higher 
density thresholds (OD-3 and OD-4; not shown in the Figure). 
Nevertheless, it appears that the thermodynamic properties of the 
gas are quite robust in AMR simulations. 

The gravitationally unstable clumps of gas identified by the 
clumpfind algorithm in the static-grid simulation at time t = 
5 Myrs are depicted in Figure 4a The corresponding results for 
the AMR runs are shown in Figures [4bf[4d| Table [T] lists the to- 
tal number and the mean size of the clumps for each simulation. 
Also listed are the fractal dimensions of the gas regions with 



4e 



number density n > 20 cm , which are plotted in Figures 

m 

For refinement by OD, the fragmentation of the CNM is 
severely underestimated. The number of clumps is roughly half 
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Fig. 4: Clump distributions (upper panel) and gas with number density n = 20 cm -3 (lower panel) for the AMR simulations compared 
to the static grid simulation. The clumps (displayed in red in the online version) were identified by a clumpfind algorithm and are 
superimposed on a grayscale density plot. 



Table 1: Number N and average size (d) of the clumps in the 
CNM; fractal dimension D of gas with number density greater 
than 20 cm -3 . 



run 


N 


(d) 


D 


static 


1025 


44 


1.45 


OD 


505 


102 


1.34 


OD-3 


481 


128 


1.33 


OD-4 


420 


157 


1.39 


CT 


1239 


32 


1.50 


RCEN 


1015 


85 


1.42 



of the number in the static-grid simulation, and the clumps are 
typically larger. The lower degree of cold gas fragmentation re- 
sults in a smaller fractal dimension (also see Figure [4f|. If the 
criteria OD-3 and OD-4 are applied, the number of clumps de- 
creases further, while their average size increases. In the case of 
criterion OD-4, a slightly higher fractal dimension is obtained, 
because the cold phase tends to fill broad, area-filling regions. 
The cooling time criterion CT yields an amount of dense clumps 
with an average size that compares well to the reference sim- 
ulation (see Figure 4c), although the degree of fragmentation 
appears to be overestimated slightly. However, we found that 
this overestimation decreases with the further evolution of the 
colliding flows and, thus, appears to be transient. Refinement 
by RCEN also reproduces the number of clumps and the frac- 
tal dimension of dense gas very well. However, there are some 
anomalously big clumps, which contribute to an average clump 
size that is systematically too lar ge. I n the plot showing gas at 
density n > 20 cm -3 (see Figure]4h|, on the other hand, such 
anomalous structures are not visible. Although refinement by 
RCEN does not overproduce gas in the cold phase (as one can 



see from the excellent agreement of the density and temperature 
pdfs in Figures 3c and [3^), there appears to be a bias toward big- 
ger clumps with this refinement method. 

In contrast to the phase separation and gas fragmentation, 
striking deviations of the turbulent flow properties in the AMR 
vs. static grid simulations become apparent. Generally, a lot of 
turbulent small-scale structure is missing in the AMR simula- 
tions. Even for the criterion RCEN, which is based on control 
variables related to turbulence, this is apparent from the con- 
tour plots of the squared vorticity modulus shown in Figure [5] 
Basically, the perturbations of the velocity field imposed at the 
inflow boundaries are quickly smoothed out in AMR simula- 
tions, so that turbulence is only produced by secondary (e. g., 
Kelvin-Helmholtz) instabilities at the collision interface in the 
central region of the computational domain. The reason is that 
all AMR cirteria, including RCEN, select relatively large fluctu- 
ations, whereas smaller perturbations are suppressed. On a static 
grid, on the other hand, the perturbations are transported from 
the boundaries to the centre and actively contribute to the pro- 
duction of turbulence. Consequently, small eddies are present in 
almost the whole domain in this case. Accordingly, the probabil- 
ity distribution of vorticity is markedly different (see Figure [6]). 
In contrast, Schm idtet al.| ( |2009| ) found very close agreement of 
the vorticity pdfs in a static-grid and an AMR simulation with 
criterion RCEN for turbulence in a periodic box with large-scale 
forcing. Our results thus indicate that the merits of different re- 
finement schemes are non-universal but rather depend on the 
properties of individual flow structures. 

4. Conclusions 

We performed two-dimensional simulations of colliding flows 
of warm atomic hydrogen with a radiative cooling function as 
source term in the energy equation. The goal of our study was 
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Fig. 5: Plots of the squared vorticity oj 2 in the static-grid simulation (a) and in the AMR simulation with refinement criterion RCEN 
(b), where the contours of the refinement levels are shown (color-coded in the online version). 
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Fig. 6: Pdfs of the vorticity modulus co for AMR-criterion RCEN 
(black curve) compared to the static grid simulation (red curve 
in the online version). 



the systematic comparison of AMR simulations, where different 
criteria for refinement were applied, to a reference simulation on 
a static grid. 

While the probability distributions of mass density and tem- 
perature are well reproduced in AMR simulations, regardless 
of the refinement technique, differences become apparent in the 
fragmentation properties of the cold gas phase. As indicators, 
we used the total number of clumps and their average size. The 
clumps were identified by a clumpfind algorithm. In addition, we 
calculated the fractal dimension of dense gas, assuming a num- 
ber density threshold of 20 cm -3 . Remarkably, the largest devi- 
ations from the clump statistics and fractal dimension extracted 
from the static-grid simulation, were encountered for refinement 
by overdensity, which is a commonly used refinement criterion 
in astrophysical AMR simulations. The deviations increase with 



the chose n density threshold. In t his regard, it is important to 
note that Hennebell e et al.| ( |2008| ) applied a density-based re- 
finement criterion, where the thresholds were chosen even higher 
than those considered in our study. Good agreement, on the other 
hand, was obtained if the cooling time or enstrophy in combina- 
tion with the rate of compression (the negative rate of change of 
the velocity divergence) were applied. 

Substantial problems with AMR became apparent with re- 
gard to turbulent flow properties. Basically, none of our AMR 
runs were able to reproduce even remotely the small-scale struc- 
ture of turbulence and the probability distributions of turbulent 
flow variables such as the vorticity modulus. This defficiency 
can be attributed to the selection effects introduced by adap- 
tive techniques. The definition of thresholds for triggering re- 
finement either selects strong local fluctuations (for example, 
large shear that gives rise to Kelvin-Helmholtz instabilities) or 
large-scale perturbations such as accumulations of mass that be- 
come Jeans-unstable in self-gravitating gas. In this respect, the 
test problem we investigated in this work is particularly tough, 
because turbulence stems from small-scale instabilities that are 
seeded by weak initial perturbations. The varying grid resolution 
in AMR simulations inevitably modulate the growth of these in- 
stabilities and, as a consequence, the production of turbulence is 
suppressed. This defficiency might be overcome by the applica- 
tion of a subgrid scale model, which transports turbulent energy 
contained in small eddies that are resolved on finer grids across 
coarser grid regions (see |Maier et al.|2009] ). 

The key point of using AMR is the computational cost for 
a given effective resolution. Indeed, Table [2] demonstrates that 
a substantial reduction of computation time is achieved with 
AMR, especially if refinement by overdensity is applied. So, 
AMR is essentially a trade-off, where fast computation has to 
be weighted carefully against the question whether the essential 
physics of the specific problem is captured. However, perform- 
ing three-dimensional AMR simulations with very high effective 
resolution is conflicting. On the one hand the reduction of com- 
putation time will definitely be even greater, on the other hand 
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there is the potential pitfall of inferring results that are properties 
of the numerics rather than the physics, since a comparison to a 
static-grid simulation is neither feasible nor desirable. 



Table 2: CPU time for the simulation runs. 



Method 


CPU time 


static 


1172h36min 


OD 


54 h 20 min 


CT 


225 h 18min 


RCEN 


290 h 57 min 
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